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abstract 


The genetic algorithm <GA) is adapted and used to obtain 

vh ctnries for methyl methacrylate 
optimal temperature histones 

polymerizations. The reaction time is minimized. while 

simultaneously requiring the attainment of design values of the 
final monomer conversion and number average chain length. The 
technique is robust, and gives near-global-optimal solutions. As 

such, it can easily be used for on-line optimizing control of free 

• T . 7 v -, a ^ Vi the reaction is 

radical polymerization reactors m which 

associated with the Trommsdorff effect. The results obtained from 
GA can be improved further if these are provided as initial 
guesses to a computer code using the Pontryagin's minimum 
principle with the first order control vector iteration method. 
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CHAPTER 1 


INTRODUCTION 

A considerable amount of research has been reported in 
the last several years on the modeling and optimization of free 

radical polymerizations exhibiting the gel or Trommsdorff 

12 '3 

effect ' . The various models have been reviewed by 0 Driscoll 

4 5 6 

and Hamielec and more recently by Achillas and Kiparissides ' 

7 

and by Mita and Horie . These models have been used in several 

g 

optimization studies which have been reviewed by Farber as well 

as by Louie and Soong . Yet, there are several unanswered 

questions. For example, Faldi et al . ^ ^ have measured the 

diffusion coefficients of methyl methacrylate (MMA) and other 

model compounds in MMA - polymethyl methacrylate (PMMA) systems 

using forced Rayleigh scattering and field-gradient NMR, and have 

inferred that the propagation rate constant, k^ (see Table 1 for 

kinetic scheme) , is not diffusion controlled, contrary to the 

general assumption used in almost all theories. They claim that 

the decrease in k with increasing monomer conversion ( x_ ) , 

p tn 

which is necessary to fit experimental data on the rates of 

polymerization, needs an explanation different from the one 

12-14 

traditionally being offered. Russell, Gilbert and coworkers 
are developing improved theories along similar lines. However, the 
earlier theories are fairly good and are still being used to 



TABLE 1 


KINETIC 


Initiation 


Propagation 

Termination by 
combination 

Termination by 
disproportionation 

Chain transfer 
to monomer 

Chain transfer 

to monomer 
via solvent 


SCHEME FOR POLYMERIZATION OF MMA 


-» 2R 


R + M =— » P„ 


P + M 1— > P , 

n n+1 


P + P ^ D (k = 0 for MMA) 

n m n+m tc 


P + P D + D 

n m n m 


P + M > P. + D 

n In 


P + S » S' + D 

n. n 

»„ fast _ _ 

S + M > S + P n 


or 

k 

P + M + D + P. 

n n 1 


2 



model, optimize and control industrial reactors, even though they 

are semi-empirical in nature. 

One group of these theories has originated from the 

molecular theory of Chiu et al . 15 Chiu et al. relate the decrease 

of the rate constants, k^ and k^_ , to the polymer concentration and 

the average molecular weight (the latter, through the initial 

concentration, [I] , of the initiator) at any time, t. Achilias 

5 6 

and Kiparissides ' related some of the parameters of this early 

model to quantities which could be measured directly using 

non-polymerizing systems. There was only a single curve- fit 

parameter, j , in their model, which was correlated to the 

initial value of the number average chain length, ^ no - The 

qualitative trends of the experimental data on the isothermal 

16 17 

polymerization of MMA in small ampoules ' (namely, sharp 

increase in x m and the weight average chain length, m w , with t, 

after the onset of the Trommsdorff effect, and the reaction 

stopping short of complete monomer conversion even though the 

reactions are irreversible, the latter being referred to as the 

glass effect) were well explained by this theory, while 

quantitative agreement was ensured by curve-fitting the value of 

the parameter. One major drawback of these versions of this group 

of theories 5, 6 ' 15 was that they could not be applied to semibatch 

reactors, or to reactors operating under non- isothermal 

conditions, and so were of little use in industry where such 

18 

operations were routinely encountered. Recently, Ray et al . and 


3 


1 9 

Seth and Gupta have presented a theory which relates the rate 

constants (as well as initiator ef f iciency^'*, f) to the current 

values of the number average chain length, 4 n - The parameters of 

this theory, 9^, 9 and 0^ (all functions of temperature, T) have 
19 

been estimated for MMA polymerization using the experimental 

data of Schulz and Harborth^ and Balke and Hamielec 17 , under 

isothermal conditions in small ampoules. The theory so 'tuned' has 

been able to explain quantitatively , experimental data on MMA 

polymerization in a 1-liter PC-interfaced, stainless steel, Parr 

reactor using idealized conditions mimicking industrial operations 

(namely, step changes in temperature and step increases in the 

2 1 

initiator and monomer concentrations ) . No additional retuning of 

the parameters was found to be necessary. This suggests that the 

theory reflects all the physico-chemical phenomena associated with 

polymerization quite well. It is to be noted that other groups of 

theories could have been modified suitably to apply to industrial 

22 

systems, but it is well recognized that almost all theories are 
about equally successful in explaining rate data and so the use of 
the relatively simple and continuous models in the group 
originating from Chiu et al . , is justified. 

18 19 

In this work, the recent theory ' for MMA 

polymerization has been used to study the optimization of a batch 
reactor. A commonly studied problem ' is to obtain the 
temperature history, T(t) (the control variable), which minimizes 
the total reaction time, t f , while simultaneously requiring the 


4 



final monomer conversion, x m £, and the final value of the number 

average chain length, to meet certain specifications (called 

'desired' values, and M nc j) • This ensures economic operation as 

well as product property requirements, and is referred to as the 

8 

minimum-time problem . A newly emerging technique, called genetic 
• 23-25 

algorithm (GA) , is used to obtain optimal solutions. This is 

an extremely robust technique and gives solutions which are quite 

close to global optima reasonably fast. Hence, this technique, 

coupled with a model which is applicable for industrial reactors, 

is well suited for use for on-line optimizing control of large 

scale MMA polymerizations (or of other similar free radical 

polymerizations), provided we can estimate the 'state' of the 

system on-line. Current experimental work along these lines is in 

progress, in which the viscosity of the reaction mass is used to 

2 6 

estimate the state of the system using inferential state 

2 7 

variable estimation . GA is then used to predict the future 
control action (temperature is being used as the control variable) 
on a supervisory- level or 'master' PC-AT, and then a 'slave' PC-AT 
implements this on a 1- liter reactor. To the best of our 
knowledge, the use of GA for predicting optimal temperature 
histories has not been reported for any system in the open 
literature on polymer reaction engineering, nor has GA been used 
for carrying out an optimally controlled polymerization 
experimentally. Most of the earlier works on the optimization of 
polymerization reactors use the far less robust Pontryagin's 


5 


27^30 

minimum principle , or the constrained pattern search 

9 

technique , to solve a variety of optimization problems as 
described in the review of Farber®, using temperature or initiator 
addition rates as control variables. These techniques are not 
suitable for use for on-line optimization work. GA, on the other 
hand, is a new and extremely powerful search technique based on 
the mechanics of natural genetics and natural selection. This 
algorithm was introduced in the mid sixties by Holland and a 
discussion of the technique and its adaptations, as well as its 
major applications are available in several books^ 4 ' ^ . It 
involves a random search over the control variable domain after 
the problem has been appropriately 'coded', usually in terms of 
'strings' or 'chromosomes' comprising of binary numbers. The best 
few solutions 'evolve' over 'generations' using techniques which 
mimic genetic evolution (hence the name) . This new technique has 
been proved to be very efficient, specially in cases where the 
objective function is flat and exhibits several local optima. The 
advantage of GA lies in the fact that it works without requiring 
much information about the system, in contrast to the traditional 
techniques which need gradients, initial guesses, etc. Hence, for 
more complex systems where the gradients cannot be easily 
evaluated and the initial guess becomes crucial, GAs lead to 
solutions which are very close to the global optima, or in fact, 
provide very good initial points to start off other techniques 
which require excellent initial guesses (e.g., Pontryagin' s 


6 


minimum principle using the first order control vector iteration 
method) . 


7 


CHAPTER 2 


FORMULATION 

Table 1 gives the kinetic scheme for MMA polymerization 
(with k^. - 0) . The mass balance equations for MMA polymerization 

in a semi -batch reactor are given by equations having the general 
form 

dx/dt = F (x,u) ; x (t=0) = x Q (1) 

where x(t) is the vector of state variables defined by 

X = [ I,M,R,S,A 0 ,A 1 , A 2 ,fZ 0 ,M 1 ,M 2 ,C m /C ml ] (2) 

and u(t) is the vector of control variables [in the present case 
it is a scalar, T(t)] 

u (t) = u (t) = T (t) (3) 

t h 

A^ and ( k = 0,1,2,...) represent the k moments of the chain 
length distributions of species P n and D n , respectively. £ m , C ml 
are additional variables to account for addition and vaporization 
of monomer after time t = 0, and are useful in the definition of 
the monomer conversion, x , for semibatch reactors. The other 
symbols are defined in the nomenclature. The exact mass balance 
and moment equations (functions, F, in Eq. 1) have been given by 

Ray et al 18 as well as by Seth and ‘ Gupta 19 (Table Bl, Appendix 
B) . These equations involve the several rate constants shown in 
Table 1, as well as the initiator efficiency, f, which quantifies 


the wastage of primary radicals, R, due to reactions not included 
in Table 1. 

The initiator efficiency, f, and the rate constants, k 

P 

and are diffusion controlled and are given by the 

1 9 

following equations 


f = f 


M 

1 + e f (T) ex P{ ?I3 w - * xef )} 


-1 


= f (x,u,p) 


(4a) 


k - 
P 


im; e P < T »(-4y) e*p{ e 13 ( 


0 - V 


ref 


>} 


-1 


= k (x, u, p) 

Jr ^ ^ 

(4b) 


"td 


T~7 — + V T) (-VT-) ex P(« - *re£) 

td, o 1 


-1 


= k td (x,u,p) 


(4c) 


with 

P 


[Sf, e p , 0 t ] 


(5) 


In this equation, is the volume of the (liquid) reaction mass, 
while ip and ^ re £ are parameters related to the weighted average 
free volume fraction as well as the molecular weight of the 
'jumping ' unit of the polymer and are functions of temperature 
(defined in Table B2) . £ I3 and ? 13 are related to the ratios of 
molecular weights of equivalent jumping units, and are constants. 
The model parameters, e f , and © t , have been tuned using the 
isothermal data of Balke and Hamielec on MMA polymerization m 
small ampoules (see Table B3 for exact expressions) . The model has 


9 


been found to be in good agreement with the experimental data on a 
® 2 0 21 

1- liter Parr reactor ' . No retuning of the values of the 

parameters, p, were found to be necessary. 

The objective function, I, used in this study is given 
by 

Min I [u(t)] . t f + VI - x mf /x md )2+ w 2 (1 ' “nf /u nd )2 (6a> 
subject to (s.t.) 


dx/dt 

= 

F(x,u) 

(6b) 

u . 

< 

u(t) r u 

(6c) 

mm 


max 

x (t) 
m 

= 

11 - M/ ^i> 

(7a) 

u (t ) 
n 

= 

<x l + U l )/( V “o’ 

(7b) 

x mf 

= 

x m (t f> 

(7c) 

^nf 

= 

W 

(7d) 


In Eq. 6, x md and n d are the desired values of monomer conversion 

and the number average chain length at t = t^, x m ^ and are the 

actual values corresponding to t = t ^ , and w and w 2 are (large) 

weightage factors. The choice of the objective function in Eq. 6 

minimizes the deviations (due to large values of w and w 2 ) of x mf 

and u r- from their desired values. The form of I used in Eq. 6 in 
nf 

which the end point requirements (constraints) are included as 

8 25 

'penalty functions', is quite popular ' . The choice, x mf = x md , 

forces the amount of unreacted monomer to be small, thus keeping 
post-reactor separation and recycling costs low. The choice, u nf 
= M nd , forces the polymer properties to be as per specifications 


10 


since several physical properties of polymers are related to the 
value of their The objective function in Eq. 6 has been used 
earlier by Sachs, Lee and Biesenberger^ 1 but with a different 
kinetic model, and by Farber and Laurence^ for styrene 
polymerization. The initial values, x q in Eq.l are given by 


[ I o ,M o ,0,0,Q,0,0,0,0,0,M o ,M o ] 


Fig.l gives the flow chart illustrating how GA, as 

applied to the present problem (Eq. 6), works. We have had to make 

24 25 

several adaptations to the conventional algorithm ' and the 

25 

computer code (SGA ) in order to solve the present problem. 

Initially (at generation number, N = 0) , a population having N 

S 3? 

chromosomes, 1.^ ; i = 1,2,...,N , is generated. Each chromosome 

chr P 

in this population comprises of a sequence of N^ a numbers (called 
substrings) which are binary representations of values of the 
control variable at equispaced points in 0 i " t fo (t fo' an 

initial estimate of t^, is to be supplied) . Each of these 
substrings, in turn, comprises of a set of N gtr binary numbers (0 
or 1) . Thus, each chromosome has N chr = N ga N str binar Y digits. The 
N chr individual binaries are generated using a random number 
generator subroutine. The binary string (sequence of N chr 

+- "h 

binaries) of the i chromosome, when decoded and interpolated 
(mapped) between the upper (u b) and lower (u ^ a) bounds of u, 
gives a digitized u-history (a set of N ga values), [u (l) ] = 


11 


1. Input Data: 

**str’ **ga* **sim' ^m* ^c' ^fo* 
S ’ “min’ “max’ 4u min’ 4 “max’ * tc ' 
"chr = V * "str 


Initialize generation 


number: 


3. Create N chromosomes, 
P 


each having binary numbers: I I Random No 


• i=12 N 

f X — X, x. f •• 0 9 « 


Generator 


4. Decode and adaptively map (Table 2) the 
chromosomes to give N digitized 

r 

u histories, each having N (decimal) 

ga 

values: 

[« (i) ] = [ u (i) (l), u (i) (2),..., u (i) (N ga ) ]; 

l = 1,2, ... , N 


Contd 


Fig. 1 Flowchart indicating the working of GA, 










Fig . 1 (contd . . . b) 


V 
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Fig . 1 (contd . . . c) 


y 
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[u ( ' (1 ) , u ( ^ (2) , . . . , u (i) (N ) ] , corresponding to that 

chromosome, as described in Table 2 (in the conventional GA 24,25 , 

k = u max an d a = u m ^ n ) . At this stage we have a set of 

chromosomes, each representing a different digitized u(t) history, 

appropriately coded in the form of a string of binaries. It 

is obvious that the minimum difference between any two digitized 

N . 

, C$ T~ -V- 

values of u is (b - a)/ (2 -1), this being the accuracy to which 

u can be determined. 

The values of u ^ (j) generated by the above procedure 
using b = u and a = u . , could fluctuate wildly between these 

two limits. This leads to significant oscillations in the optimal 
u-histories, and is both undesirable and non -implement able . In 
order to reduce these oscillations, further constraints are 
clamped on to the values of u^ (j) in the present study, so that 
neighbouring values of u do not differ by more than some 
prescribed values, Au m ^ n and Au max - Thus, 


or 


iu min s Au<l) ,j) + s iu max 


u<1) ( 3> + iu min s u<1 ’ (j+1) 3 u<1> (j) + iU max 


[9) 


( 10 ) 


(i) 


where Au . is a negative number. Thus, the first value, u (1) , 

corresponding to t = 0, is determined randomly to lie between u max 

and u . , while all subsequent values are determined randomly 

min' 

within a smaller range around the previous value. This procedure 
is being called adaptive mapping. The accuracy (minimum 
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TABLE 2 


DECODING AND ADAPTIVE MAPPING PROCEDURE FOR 


N =2 
ga 


N chr ‘ % N str ' 2N str • a s “ 3 b 


Example : 


= t 10011. . .0 ; 


11100 ... 0 ] ; i = 1, 2, . . . ,N 


N str binaries N str binaries 


'Decode' each of the (two sets of) binary numbers into decimal 
numbers, d.^ and , using, for example, 

N . -1 N , -2 0 

, , _ str _ 0 str _ „ 

d 1 =1x2 +0x2 +...+0x2 

Now, using the 'mapping' obtain the digitized u-history 


[ u (i) (j) ] = [ u (l) (1) , u (l) (2) ] 


a + d . x VAL 
1 


; i = 1 , 2 ...... N 

r 

j = 1,2 


where 


VAL = (b - a) / (2 b - 1) 


* a = u min' b " u max ; for j 1 

a = u (i) ( j - 1) + Au min , b = u (l) (j - 1) + Au max 

for j > 1; s.t. u min * a, b s u max 
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difference between values of u) of internal points is observed to 

be higher than for the first (t = 0) point, if lAu - Au . I < 

c 1 max mm 1 

I u - u . I . 

1 max mm 1 

The decoded and adaptively mapped, discretized values of 

u are curve- fitted piece-wise (splines) to obtain a continuous 

function, u^^ (t) . A piece-wise cubic Hermite subroutine is used 

to do this. This continuous function is again digitized to give 

N .( - N ) values of the control variable, [U^ (j); j = 
sun ga 

1,2,...,N . ]. The generation of several additional, intermediate, 
s im 

discretized values of u^ is necessary for integrating the model 

differential equations (Eq. 1 and Tables B1-B3) . 

The digitized temperature history, [U ^ (j); j 

t h 

1,2,...,N . ], corresponding to the l member of the population, 
sim 

3 3 

is used in a Gear subroutine (D02EJF in the NAG library) to 
integrate the balance equations, starting with the initial 
conditions in Eq. 8 and continuing till t = t^ . The program 
stores the values of each of the state variables, x (i) (j), at 
every intermediate value of t, such that there are sets of x. 

The value of I ^ at each of these storage locations is computed, 


and the location, tj3^, of the minimum of l' ,x ' as well as the 
minimum value itself, 1^, are obtained by search. Evidently, t fQ 
should be chosen large enough so that I m ^ n occurs in 0 ^ t t^ Q 
for all i. The integration of the balance equations and the 
location of 1^^ for each of the N p chromosomes is carried out. 


(i) 
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One additional point needs to be emphasized. The 
2 5 

computer code, SGA , which has been used in this study after 
modification maximizes a 'fitness 7 function, F^^ , rather than 
minimizes an objective function, I. Hence, we define a fitness 
function 


F (l) si/ ( 1 + I ( ! } ) 

mm 

(i) 


( 11 ) 


and maximize its value (wherever I v ' is to be minimized). 

The next step in GA is to have 'reproduction' in the 

population of chromosomes. A mating pool is first formed. In this 

pool, priority is given to those chromosomes which have higher 

fitness values. The essential idea is to pick out the 

above-average strings in the current population and include 

(multiple) copies of these in the mating pool in a probabilistic 

manner. It is here that the principle of natural selection 

(survival of the fittest) comes in action. The principle of 

proportionate reproduction is used. The probability of selecting 

, , . . N , . , 

the l chromosome in the mating pool is F 1 / Yp IF ' 1 ■ A 

i=l 

'roulette-wheel' (whose circumference is marked for each 

chromosome proportionate to its fitness value) is spun N times. 

ir 

In each spin, the chromosome corresponding to the location of the 

roulette -wheel pointer is copied into the pool. This thought 

24 25 

experiment is implemented using random numbers ' 

After the mating pool is created, crossover and 

mutations take place to produce . the new population (next 

generation) . These operations take place at the chromosome 
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(binary) level . Two chromosomes are selected randomly from the 
mating pool, a crossing site is selected (randomly again), and 
portions of the chromosomes before and after the crossing site are 
exchanged. For example, for seven-bit chromosomes with crossing 
site after the third binary, the crossover is described by: 


100 | 1111 

110 I 0100 
(old generation) 


100 0100 
110 1111 

(new generation) (12) 


While performing crossovers, only NpP c chromosomes are crossed, 
the remaining being left untouched (p is referred to as the 
crossover probability) . 

Another operation, called mutation, is also used to 

improve the next generation. The mutation operator changes a 

binary number from 1 to 0 or vice versa, with a probability, p . 

This operation is carried out for each of the NpN c ^ r bits in the 

24 25 

population, again using appropriate random numbers ' . The need 

for mutation leads to a local search around the current solution, 

25 

and helps maintain the diversity of the population 

The random crossover procedure discussed above leads to 
a preponderance of crossovers in the (inactive) range, tju^ =s t s 
t f , if the guess value of t fQ supplied to the computer code is 
too large. This procedure, thus, needs to be adapted so that 
crossovers take place in a t-domain (horizon) which becomes 
smaller over generations. What is done is to limit crossovers to 
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is the best (minimum) 


0 < t s t 

mm, best 
of the Np values 

factor supplied 

experimentation) . 

S) is given by 


x (1 + S) , where t . , 

mm, best 

°f t t riin in any generation, and S is a safety 

to the program (obtained by numerical 

The string length corresponding to t . , (1 + 

min^Dest 


N chr N chr t min,best^ 1 + / t f Q 


(13) 


where N c ^ r is an (next higher) integer. The region in which 
crossovers take place would decrease from generation to generation 
as ^min best decreases. Such an adaptation of the conventional GA 
can be used to advantage for any minimum time optimization 
problem, and provides an automatically narrowing crossover 
horizon . 


The optimal solutions generated using GA can be compared 

with those obtained (for the same objective function, constraints 

27-30 

and model equations) from Pontryagin's minimum principle 

using the first order control vector iteration technique (referred 

to as PI) . The algorithm to be used is summarized in Table 3, and 

34 

is an adaptation of that used by Vaid and Gupta and Ray and 
3 5 

Gupta earlier. 
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TABLE 3 


FORMULATION OF THE OPTIMAL CONTROL PROBLEM USING 
PONTRYAGIN'S MINIMUM PRINCIPLE 27 "" 30 WITH FIRST ORDER 
CONTROL VECTOR ITERATION METHOD 
Optimization problem 

Max I [u (t) ] = G [x(t f ) ] 
s . t . 


Here , 

G [x(t f )j = 


dx/dt 


u . 
mm 


= F (x, u) 

U (t) 2= U 


max 


t f + w 1 (1 - 


x £ /x ,) 
mf' md 


+ w 2 (1 


^nf^nd 


)* ] 


Procedure 

1) Guess u(t) = (t) ; 0 s t £ 

2) With this u(t), integrate the state variable equations 

dx/dt = F(x,u) to obtain x(t) ; 0 £ t ^ t^, 

with t f obtained by solving 

H(t f ) = (3G/3x) F | t = = 0 

3) With the values of x(t) and u(t), integrate the 
adjoint equations backwards from t = t^ to t = 0, 

dk T /dt = - (SH/3x) ; A T (t f ) = 3G/3x| t = t 

where H = X^F 

4) Correct u(t) by, 

u(t) = u(t) + e (3H/3u), c > 0 

5) Perform a single variable search (on s) to 
generate several u(t) , obtain I for each of these 
histories (after integrating the state variable equations 
for each case), then obtain c Qpt corresponding to 

Contd. . . 3b 


, } Ti. if? 

i S. T„ KANtbtt 

, millllll llllw I'lll'linrr 

0 4 O C 
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Table 3 (contd. . .b) 

the maximum value of I (note that SH/du is not updated 
during this search) . 

6) Update u(t) by, 

u neW (t) = u old (t) + e opt (SH/Su) 

and return to step 2 . 

7) Iterate until convergence is attained. 
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CHAPTER 3 


RESULTS AND DISCUSSION 

Several checks were made to ensure that the computer 

code prepared was free of errors. The code was used to generate 

the monomer conversion and the number average chain length for 

different isothermal conditions. These were found to match the 

1 9 

results of Seth and Gupta and are shown in Fig. 2 for an initial 
initiator concentration, [I] , of 25.8 mol/m 3 (reference value). 

These results were generated with a value of TOL of 10 -7 in the 
code D02EJF, and no significant differences were found upon 
decreasing the value of this parameter. This check indicated that 
the simulation part of our code was free of errors, and also 
provided results which could be used to explain optimal histories 
qualitatively. The next check was on the correctness of the 
optimization part of our program. From Fig. 2, it is clear that if 
we use 

x , = 0.0134 
md 

JU . = 2365 (14) 

nd 

6 0°C £ T(t) s 90°C 

the optimal temperature history would be isothermal at 60°C (any 
higher temperature would give lower values of h n £ while 
simultaneously giving higher x m f) • Also, under these isothermal 
conditions, the value of t^ would be 230.77 s. Similarly, for 




0.4926 


X md " 

M nd = 532.16 (15) 

6 0°C s T (t ) £ 90°C 

the optimal T(t) would be isothermal at 90°C (see Fig. 2), with t f 
= 969.67 s (any lower temperature would lead to higher values of 
d nf , while simultaneously leading to lower x mf ) . The optimization 
problems described in Eqs . 14 and 15 were solved using the GA 

computer code (using the parameters of Table 4) and in both cases, 
the expected optimal temperature histories were obtained (in 8 
generations for Eq. 14 and 3 generations for Eq. 15, shown in 
Figs. A1-A2 in Appendix A) . A similar check was made for the 
computer code using Pontryagin's principle with the first order 
control vector iteration method (PI) . The starting guess for this 
technique was T^°^= 90°C (for Eq. 14) and T^°^= 60°C (for Eq. 15). 
Again, the expected isothermal optimal histories were obtained 
(Figs. A1-A2 in Appendix A) in 2 and 8 iterations (for Eqs. 14 and 
15, respectively) . These checks gave confidence on both our 
computer codes, GA and PI. 

The optimization program using GA was now run for 


U , = 1800 

^nd 

34 

These values are quite close to those used by Vaid and Gupta, as 
well as other workers, and are being used as reference values to 
illustrate the working of GA. Fig. 3 shows how the optimal 
temperature history (the best for each generation) evolves over 


25 



500 10' 

Pi g . 3 : Evolution of temperature histories towards the optimal one , wt 
generation number, N g , corresponding to x md = 0.94, ^ 

= 1800 (for parameters of Table 5 ). Arrows indicate the end 
points of corresponding curves . 





2500 


o T (t) for the GA15 run , as well as 
opt 





TABLE 4 


PARAMETERS USED FOR REFERENCE RUN 

GA Parameters 


N = 100 
P 

N str= 7 

V 10 

N . = 100 

sim 

^ u min ' U max ^ 

[ Au . , Au ] 

mm max 


[ 60 , 90 ] ; °C 

= [ -15 , +15 ] 


o 


C 


P c = 0.99 

p = 0.000009 
*m 

S = 0.2 

t c = 4 0 00 s 
ro 

w = w 2 = 2.5X10 5 

In addition, the value of the parameter, RS, 
binaries is 0.9. 

Design Parameters 


,25 

used 


for generating 


x md ~ 
^nd = 


0 . 94 
1800 

25.8 mol/m 3 
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generations. Very little improvement takes place after about 12 

generations and so, results for N > 12 are not shown. The CPU 

a 

time for generating these results was 15.8 s on a DEC 3000 axp. 
The variation of and x m with time, using the optimal 

temperature history (for N g = 12) shown as GA15 (15 indicating 

^ u min anc ^ ^ u max an< ^ +15 ° c ) in Fig. 3, are shown in Figs. 4 

and 5 (by solid lines marked GA15) . Some amount of oscillations 
are observed in the optimal temperature history (Fig. 3) , which 
could be reduced by changing some of the parameters in Table 4 
(see later) . Fig. 6 compares the optimal history (curve GA15; same 

as for N = 12 in Fig. 3) with that obtained using Pontryagin's 

3 

minimum principle (curve PI, obtained by starting with T^°^ = 90°C 

and converging in about 8 iterations) . The values of the objective 

function, I, for the GA15 and the PI cases are found to be 2008.36 

and 2016.96, respectively. The two histories are also observed to 

24 25 27-30 

be fairly close to each other. It may be noted ' ' that 

both GA and PI lead to near-optimal solutions only and become 
very sluggish as the optimal history is approached. The agreement 
in Fig. 6 is, thus, extremely good (perhaps fortuitously so) . It 
is interesting to observe from Figs. 3 and 4 that optimal 
operation requires relatively low temperatures [leading to 
relatively high values of fi n (see Fig. 4)] , followed by a gradual 
increase in T(t) [associated with some fall in #i n ] to its maximum 
value of 90°C . The value of u n builds up to its desired value by 
exploiting the gel effect near the end, this being exhibited as a 
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Fig. 6 : T^(t) corresponding to the conditions of Fig. 3 using the PI ( dotted ) 
and GA ( solid ) techniques . GA15 corresponds to the reference run 
( Table 5 ) while GA30 corresponds to [ AT . , AT ] = + 30°C 

v 7 1 mm max 

( all other parameters same as given in Table 5 ) . Arrow indicates t for PI . 


31 




sharp increase in n n (t) and x^t) near t « t £ . The sudden increase 

^n^ t0 ; '‘ ts value of is a characteristic of almost 
all optimal solutions obtained in our study, and emphasizes the 
need for model-based on-line optimizing control in the period 
prior to the onset of the gel effect. 

Computations were carried out using Pontryagin's minimum 
principle (first order) for the conditions described in Eq. 16, 
but using the initial guess, T (o) (t), different than that used for 
generating Figs. 4-6 [i.e., isothermal T (o) (t) different from 
90°C, see Figs. A3-A5 in Appendix A]. It was found that the (sub-) 
optimal temperature histories were very sensitive to the initial 
guess, and that there was only a very narrow window of the initial 
guess for which converged solutions were obtained which were 
similar to GA15. A similar acute sensitivity to the initial guess 
history was also observed by Vaid and Gupta^ 4 who used a similar 
algorithm but solved a slightly different optimization problem. In 
fact, we obtained different (sub-) optimal temperature histories 
on using different (t) outside of the narrow window. In each 
case, the and (i n ^ were very close to their desired values, 
while the values of I differed slightly. This could be because ,of 
two possible reasons. Firstly, the value of I t is relatively 
insensitive to T ^.(t), and secondly, there could be several 
shallow, local minima, and PI converges to (near) these, depending 
on the initial guess, (t) , provided. Which of these two causes 
leads to the ineffectiveness of the PI technique is not clear, nor 
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to reach 


is this answer too important. However, GA is known 24 ' 25 
the global optimum, and is robust, and so we believe that its 
solution is the true' one. The PI technique also converges to the 
solution provided by GA if we start from T (o) = 90°C, or use 

T (t) somewhat similar (but not identical) to the optimal 
history provided by GA. 

This drawback of the PI technique can be overcome by the 
use of GA to first generate near-optimal solutions which can be 
provided as initial guesses, to be improved upon by using PI. We 
believe that GA followed by Pi is a superior combination than the 
first order Pontryagin technique followed by the second order (P2) 

9 Y „ 3 q 

technique " ' , in which second order derivatives are required. 

Fig. 7 shows the improvement of the optimal T(t) using the GA15-P1 
combination. The value of I of 2008.71 corresponding to GA15 is 
reduced to 3980.34 using the PI technique. Similar improvements in 
the value of I have been found in other cases of GA + PI tried in 
this study (detailed results can be provided on request) . 

We now study the effect of varying the parameters (Table 
4) used in GA. Details of the parameters which are varied one at a 
time, keeping all others at their reference values (Table 4), are 
given in Table 5. Fig. 8 shows that in the initial region (low t) , 
the optimal temperature history is somewhat sensitive to the 

number, N , of chromosomes in the population (curves 1 and 2) . 

P 

However, changing N could lead to oscillatory behavior in 

Jr 

T „ (t) , which needs to be dampened by changing some other 

opt ' 
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TABLE 5 

SOME DETAILS CORRESPONDING TO T Qpt ( t ) SHOWN IN FIGS. 8-1 



For achieving convergence 






parameter (e.g., Au min , Au max ) simultaneously. In fact, the 
reference values of the parameters (run GA15) have been chosen 
such that the oscillations are minimized for this run. Similar 
oscillatory behavior is observed (curve 3, Fig. 8) in the initial 
region by increasing from 7 (ref) to 14. It is clear that any 

change made to improve the accuracy of results leads to more 
oscillations in the initial region, and its effects need to be 
dampened out. A similar conclusion is obtained on studying curves 


4 and 5 in Fiq. 9 and curve 8 in Fig. 10. As N 

ga 

Top, ( t ) oscillates considerably, to the extent that 


is increased, 

I . worsens . 
opt 


The effect of increasing the mutation probability is similar. 


Decreasing the crossover probability from 0.99 to 0.98 (curve 7, 


Fig. 10) does not lead to oscillations, but worsens I t slightly. 
The effects of decreasing and changing RS are also shown 

(curves 8 and 9) in Fig. 10. 


Fig. 11 shows the effect of varying the parameter 

characterizing one of the - adaptations of the conventional GA, 

namely, use of Au min and Au max as constraints. These were 

introduced to dampen oscillations in T^^t), as well as to ensure 

implementability of the optimal history in industrial systems. The 

actual values of AT . and AT to be used should really be 

min iiia. a 

decided by the heat transfer limitations of the reactor, but these 
have been considered as parameters and chosen somewhat arbitrarily 
here to study their effect. It is observed that increasing the 
range of AT from ± 15 to ± 30 °C leads, as expected, to more 
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] and S on the optimal temperature history . 

1 = + 20° C and + 30° C , respectively 





i 


oscillations, and to a worsening of I It is interesting to 
compare curve 11 (for -30 ^ AT 2 30°C, i.e., no constraint is 
operative on the temperature of a neighboring point, except 60°C s 
T ^ 90°C) with the results from the PI technique [with (t) = 
90°C ] where no constraint on the temperatures of neighbouring 
points are operative. Fig. 6 shows that curve 11 (renamed GA3 0) 
does not compare as well with curve PI quantitatively , as does the 
GA15 results due to the oscillations present in GA30. Use of a 
damping mechanism through AT m j_ n and AT max' thus, appears 
justified . 

Fig. 11 also shows (curve 12) the effect of increasing 
the safety factor, S, a parameter reflecting another adaptation we 
have made in the conventional GA. Increasing S leads to a larger 
domain in which crossovers are permitted, and slows down the rate 
of convergence (note that GA, too, becomes sluggish as the optimal 
history is attained, and the results in Figs. 8-11 are all 
near-optimal in that sense) . 

The general conclusion from this parametric sensitivity 
study is that we need to experiment with the several parameters to 
obtain good, sub-optimal u-histories with GA. Since the histories 
are global (sub-) optimal solutions, we can follow up GA with the 
first order Pontryagin (PI) technique to get good final results. 
This combination exploits the best features of both these 
techniques . 

Figs. 12 and 13 show the effect of varying the 'design' 
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t , s 


Fig. 14 : |i (t) corresponding to the optimal temperature histories given in Fig. 13 . 
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variables, x md , h nd and [I] Q . Fig. 12 (curve 13) shows that 
somewhat lower initial temperatures and slower rates of rise of 
T(t) are required to obtain higher final values of the monomer 
conversion if we wish to keep unchanged. The presence of 

oscillations in (t) indicates that the reference values of the 

parameters used are not appropriate to generate the results for 
this case, and need to be 'retuned' if we wish to have better 
results. Fig. 13 (curves 14 and 15) shows how the increase of 
T t (t) should be delayed to give higher ^ nf products. Fig. 14 
shows the delayed gel effect helping achieve higher h n £ products. 
The effect of decreasing the initiator loading, [I] , is also 
shown in Fig. 13 (curve 16) . Higher temperatures are necessary 
with lower [I] Q to speed up the reaction, so that t f is minimized. 

The general trends observed in all these cases is that 

optimal temperature histories for MMA polymerization are such 

that, initially, we have almost constant • This is followed by a 

period during which n decreases (as T goes up) . Finally, the 

gel-effect occurs which leads to a relatively rapid increase in 

x and ijl to their desired values. The temperatures in the pre-gel 
m n 

effect region are quite important, particularly since rapid 


changes in T after the onset of the gel effect are not easy to 
implement. This points out the need for using model-based on-line 
optimizing control. It is difficult to predict the qualitative 


trends ot T 0 p t (t) intuitively 
in Fig. 2, and this emphasizes 


using the isothermal results shown 
the importance of such quantitative 
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studies . 


The variation of the polydispersity index (PDI, see 
Nomenclature) of the polymer with time, under optimal conditions 
(GA1 5) , is shown in Fig. 15. The final value of the PDI is 
observed to be substantially lower than that of polymer produced 
under isothermal conditions, which is a blessing in disguise, 
since reduction of the PDI was not envisaged in our optimization 
problem (Eq.6) . There appears to be some controversy in the 
literature regarding whether the minimum time problem ensures 
simultaneously , minimum PDI . Our results indicate substantial 
lowering of the PDI. 
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CHAPTER 4 


CONCLUSIONS 

A robust optimization technique, Genetic algorithm, has 
been used in this study to obtain global optimal temperature 
histories for MMA polymerization. These can be improved further by 
using the first order Pontryagin method. The technique can easily 
be used for on-line optimizing control of experimental reactors. 
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CHAPTER 5 


SUGGESTIONS FOR FUTURE WORK 

The rate of liquid monomer and initiator addition 
Rli (t) was not considered as a control variable in the 
present problem. This may be included in future as a 
control variable along with T(t) to cut down the 
reaction time further. 

The algorithm with proper tuning of its computational 
parameters can be used for on-line optimizing control of 
MMA polymerization process (work along these lines is in 
progress) . 
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t , S 


Fig. A1 : Evolution of temperature histories towards the optimal corresponding to Eq (14) 
(a) using GA , generation to generation, (b) using PI , iteration to iteration . 
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Fig. A2 : Evolution of temperature histories towards the optimal corresponding to Eq (15) 
(a) using GA , generation to generation, (b) using PI , iteration to iteration . 
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TABLE B1 


model equations for mma polymerization in semibatch reactors 


19 
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Contd. . - b 
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Table B1 (Contd. . . b) 
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TABLE B2 


CAGE, GEL AND GLASS EFFECT EQUATIONS 19 
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Table B2 (Contd. . . b) 
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TABLE B3 


PARAMETERS USED FOR POLYMERIZATION OF 

19 

MM A 
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d 
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p, O 


“td, o" 


tc 


k. 
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E 


td 


(MW ) = 
m 


(MW ) = 
s 


966.5 - 1.1 (T - 273.1) kg/m 3 

1200 kg/m 3 

844.18 - 1.07165 (T - 323.1) kg/m 3 (Benzene) 

0.58; for AIBN 

1.053X10 15 s" 1 ; for AIBN 
4.917xl0 2 m 3 /mol-s 

9.8xl0 4 m 3 /mol-s 

0.0 

0 . 0 
k 

P 

0 . 0 

128.45 kJ/mol ; for AIBN 
18.22 kJ/mol 
2.937 kJ/mol 
0.10013 kg/mol 
0.07811 kg/mol 

Contd. . . b 
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Table B3 ( Contd. . . b) 


(MWj ) = 0.06800 kg/mol; for radicals from AIBN 


Constitutive Parameters for the Cage, Gel and Glass Effects 
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T< (105+273 .1) K 

Correlations used for Q c , 6 

r p 

e 19 
' 0 t 


log Q [10 3 S f (T) , m 3 mol 1 ] = 2.016 x 10 2 - 1.455 x 10 5 (l/T) 

+ 2.70 x 10 7 (l/T 2 ) 

log [e (T) , s] = 8.03 x 10 1 - 7.50 x 10 4 (l/T) 

XU |p 

+ 1.765 X 10 7 (l/T 2 ) 

log i 0 [©t ( T ) ' = 1-241 x 10 2 - 1.0314 x 10 5 (l/T) 

+ 2.2735 x 10 7 (l/T 2 ) 
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